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Massive spin-2 fields on black hole spacetimes: 
instability of the Schwarzschild and Kerr solutions and bounds on graviton mass 

Richard Brito/'H Vitor Cardoso, 1 - 2 ' 3 -Q and Paolo Pani x ' 4 >[l 



1 CENTRA, Departamento de Fisica, 
Universidade Tecnica de Lisboa - UTL, Avenida 



Instituto Superior Tecnico, 

Rovisco Pais 1, 1049 Lisboa, Portugal 



Perimeter Institute for Theoretical Physics Waterloo, Ontario N2J 2W9, Canada 

3 Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA. 

^Institute for Theory & Computation, Harvard- Smithsonian CfA, 60 Garden Street, Cambridge, MA, USA 

Massive bosonic fields of arbitrary spin are predicted by general extensions of the Standard Model. 
It has been recently shown that there exists a family of bimetric theories of gravity - including mas- 
sive gravity - which are free of Boulware-Deser ghosts at the nonlinear level. This opens up the 
possibility to describe consistently the dynamics of massive spin-2 particles in a gravitational field. 
Within this context, we develop the study of massive spin-2 fluctuations - including massive gravi- 
tons - around Schwarzschild and slowly-rotating Kerr black holes. Our work has two important 
outcomes. First, we show that the Schwarzschild geometry is linearly unstable for small tensor 
masses, against a spherically symmetric mode. Second, we provide solid evidence that the Kerr 
geometry is also generically unstable, both against the spherical mode and against long-lived super- 
radiant modes. In the absence of nonlinear effects, the observation of spinning black holes bounds 
the graviton mass fi to be fi < 5 x 10~ 23 eV. 



I. INTRODUCTION 

The feebleness with which exotic particles- such as 
those predicted in several extensions of the Standard 
Model [l|-y] or in modified theories of gravity [j- cou- 
ple to ordinary matter, lies at the heart of the difficulty 
to detect them. Extra fundamental fields may couple to 
Standard Model particles in various ways, which makes 
it challenging to exclude, or possibly detect new effects. 

Fortunately, the equivalence principle guarantees that 
all forms of matter gravitate. Therefore, it is no sur- 
prise that extra fundamental fields - specially if extremely 
light by Standard Model standards - can strongly affect 
the dynamics of selfgravitating compact objects, such as 
black holes (BHs) and neutron stars. The equivalence 
principle, together with the fact that BHs are vacuum 
solutions, guarantees that all forms of matter, including 
exotic matter, interact with BHs in the same universal 
way. One is thus offered the intriguing possibility of using 
the growing wealth of observations in high-energy astro- 
physics J5rt7[ to put physics beyond the Standard Model 
to the test. 

There is a vast literature -which we will not attempt 
at summarizing- on the gravitational interaction of fun- 
damental scalar fields Q. Of more direct interest to us 
are recent efforts to use BHs as particle-physics labora- 
tories, through which one can constrain the mass of the 
QCD axion, of strin gy p seudoscalars populating the so- 
called axiverse [l|, l9J4Tl| and the hidden C/(l) sector of 
the Standard Model H H ES EH- In addition to their 
phenomenological relevance, such studies have revealed 



unexpected aspects related to the dynamics of these fields 
in curved spacetime. 

In this paper, we take a further step in this enterprise 
by investigating the dynamics of massive spin-2 fields 
propagating on a BH spacetime. 



Executive summary 

For the reader's convenience, we summarize here the 
structure of the paper and our main results. To put 

Hi is devoted to a generic 
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our work into context, Section 
discussion on massive gravity 

ories |15l - [l7| . on the dynamics of spin-2 fields on curved 
spacetimes and their possible imprint in gravitational- 
wave and BH physics. We also discuss how ultralight 
spin-2 fields are expected to trigger strong superradiant 
instabilities (l8l— [2lj| in massive BHs. 

In Section IIIII we review Fierz-Pauli theory [22| for a 
linearized massive spin-2 field pro pag ating on flat and 
curved backgrounds (see also Ref. [23|). The linearized 
field equations for a massive spin-2 fluctuation propagat- 
ing on curved spacetimes are given in Eqs. (|26p - (|28p and 
we also show how they can be consistently obtained in 
bimetric and massive gravity. 

Within this context, Sections IIVI and [V] are devoted 
to a complete analysis of the linear dynamics on a 
Schwarzschild BH. In Sec. IIVI we focus on the monopole 
mode that corresponds to the scalar polarization of a 
massive graviton. We find a strongly unstable, spher- 
ically symmetric mode, which was also discussed very 
recently in Ref. [2J| • Thus Schwarzschild BHs are unsta- 
ble in these theories and we show that the inclusion of 
a cosmological constant makes the Schwarzschild-de Sit- 
ter BHs even more unstable. Furthermore, in Sec. |V] we 
derive the full master equations for the axial and polar 
sectors. We find that the spectrum supports quasinor- 



mal modes (QNMs) and quasibound, long-lived states for 
any non spherically symmetric mode and we compute the 
spectrum numerically. 

In Sec. lVII we extend our analysis to stationary and ax- 
isymmetric BHs, namely to the Kerr metric. In general, 
the radial and angular part of the perturbation equa- 
tions on a spinning geometry are challenging - if pos- 
sible at all - to separate within the standard Teukol- 
sky approach [2|| [26[ . The same obstacle is encoun- 
tered for massive spin-1 (Proca) perturbations of a Kerr 
BH. The problem has been recently solved within a 
slow-rotation framework j27H29l | in the frequency do- 
main [10|, [H| and also usin g fu ll-fledged numerical evo- 
lutions in the time domain [30] - We have extended the 
technique of Refs. piJ, LLL( to the case of massive spin- 
2 perturbations [see also [31( for the case of gravito- 
electromagnctic perturbations of Kerr- Newman BHs] . 

We derive the perturbations equations to first order 
in the BH angular momentum. In principle, this proce- 
dure can be extended to any order. To first order, the 
eigenvalues of the system are described by two indepen- 
dent sets of equations (one for each parity) and for each 
harmonic index. By solving the first-order equations, we 
have found strong evidence for the existence of unstable 
modes in the spectrum. This instability is different from 
that affecting Schwarzschild BHs and it is associated to 
nonspherical modes which becomes unstable above a cer- 
tain BH angular momentum. The instability can be four 
orders of magnitude stronger than in the Proca case and 
up to seven orders stronger than in the massive scalar 
case. Our results provide strong indications that massive 
spin-2 fields trigger the strongest superradiant instability 
in vacuum BH solutions. 

Although a second-order analysis would be necessary 
to describe superradiancc consistently, a first-order ap- 
proximation is generally sufficient to give accurate re- 
sults well beyond its regime of validity [101 ] . Including 
second-order effects would be an important - and tech- 
nically challenging - extension of our work. The unsta- 
ble, spherically-symmetric mode active for Schwarzschild 
BHs is unaffected by rotation, at first order. Thus, we 
present two mechanisms by which Kerr BHs are rendered 
unstable in massive theories of gravity. 

Several technicalities are discussed in the Appendices 
and in publicly available Mathematica notebooks [32| . 
In Appendix [D] we generalize Detweiler's calculation of 
the unstable massive scalar modes of a Kerr BH |33[ to 
the dipolar axial sector of massive spin-2 fields to first 
order in the BH angular momentum. 

We conclude in Sec. I VIII with some phcnomenological 
implications and with possible future extensions of our 
results. 



II. MASSIVE SPIN-2 FIELDS AND STRONG 
GRAVITY 



A. Massive gravitons? 

Higher-spin fields are predicted to arise in several con- 
texts [3J-[3a]. The motivation to investigate their gravi- 
tational dynamics is twofold. The first reason is concep- 
tual and is tied to a renewed interest in massive gravity 
and bimetric theories of gravity. It is known since the 
work of Fierz and Pauli that at the linear level there is 
only one ghost- and tachyon-free, Lorentz-invariant mass 
term that describes the five polarizations of a massive 
spin-2 field on a flat background [22|. However, in the 
zero-mass limit the Fierz-Pauli theory does not recover 
linear general relativity due to the existence of extra de- 
grees of freedom introduced by the graviton mass. In 
the massless limit the helicity-0 state maintains a finite 
coupling to the trace of the source stress-energy tensor, 
modifying the Newtonian potential and hence yielding 
predictions which differ from the massless graviton the- 
ory [2j, [37H39J . rendering the theory inconsistent with 
observations. This is known as the vDVZ discontinu- 
ityHEI. 

To overcome this difficulty Vainshtein [42[ argued that 
the discontinuity present in the Fierz-Pauli theory is an 
artifact of the linear theory, and that the full nonlin- 
ear theory has a smooth limit for m g = hfi — > 0. He 
found that around any massive object of mass M, there 
is a new length scale known as the Vainshtein radius, 

Ty ~ (M/(m,gMp)) . The nonlinearities begin to dom- 
inate at r < ry invalidating the predictions made by the 
linear theory. This is due to the fact that at high ener- 
gies the helicity-0 mode of the graviton, responsible for 
the discontinuity, is strongly coupled to itself and be- 
comes weakly coupled to external sources. However, it 
was believed until recently that Lorentz-invariant non- 
linear massive gravity theories were doomed to fail due 
to the (re) appearance of a ghost-like sixth degree of free- 
dom [43- This was studied by Boulware and Deser who 
showed that in nontrivial backgrounds there are 6 degrees 
of freedom, where the extra degree of freedom was shown 
to be a ghost scalar, known as the Boulwarc-Dcser ghost. 

More recently, a two-parameter family of nonlinear 
generalizations of the linear Fierz-Pauli th eory was pro- 
posed by de Rham, Gabadadze and Tolley (l3 - fl4j and it 
is usually referred to as "nonlinear massive gravity" [see 
Ref. [23[ for a review]. When linearized on a flat back- 
ground, nonlinear massive gravity has so far proved to 
be ghost-free order by order (but see Refs. [44[ for recent 
counterarguments and Ref. [45j for some tight constraints 
on the theory in the decoupling limit). The extension of 
the theory to generic nonflat backgrounds appears to be 
also ghost-free [4y, |47[. On the other hand, it has been 
recently shown that the very same combination that re- 
moves the Boulwarc-Deser ghost is also responsible for 
the existence of superluminal shock- wave solutions which 
render the theory acausal [48j |. 

Furthermore, the healthy interaction term that pre- 
vents the theory to propagate ghosts has been also gen- 
eralized to bimetric theories of gravity, i.e. to theo- 



ries which propagate two dynamical spin-2 fields [15I - TT7T ] . 
These theories can also describe a massive spin-2 field 
coupled to standard Einstein gravity [49| and they re- 
duce to nonlinear massive gravity when one of the fields 
is nondynamical |50|. 



B. Gravitational- wave searches and astrophysics 

The second motivation to investigate massive spin-2 
fields is of a more practical and phcnomcnological nature. 
Advanced gravitational-wave detectors will begin opera- 
tion in a couple of years and the first direct detection of 
a graviton on Earth is expected to take place within the 
next decade. Current constraints on the graviton mass 
from pulsar observations already provide compelling ev- 
idence that gravitational waves are indeed emitted when 
two objects merge [51|. A hypothetical massive grav iton 
would affect the decay rate of the orbiting pulsar |52|, [53[ . 
The Hulsc- Taylor pulsar provides a stringent limit on the 
mass of the graviton @,/i<7.6x l(T 20 cVE]. 

However, even with these tight constraints in place, 
the Yukawa-like potential of a hypothetical graviton 
mass would be responsible for a deformation of the 
gravitational-wave signal during its journey from the 
source to the observer. In other words, a small graviton 
mass may not affect the inspiral of a binary system to a 
significant extent (including the changes in period of bi- 
nary pulsar), but introduces nontrivial dispersion which 
acts over several Compton wavelengths, ~ /z _1 . This 
peculiar effect can leave a signature in the gravitational 
waveform. Because any putative gravitational-wave de- 
tections will occur with very low signal-to-noise-ratio, an 
accurate knowledge of these effects may be important, in 
the sense that accurate templates are required to detect 
extra polarizations without introducing bias |55l l56j [see 
Ref. [57| for a recent review] . 

In summary, gravitational waveforms for inspiralling 
objects emitting massive gravitons are necessary. There 
are several ways to deal with this problem, e.g., full non- 
linear simulation, slow-motion expansions or perturba- 
tive expansions around some background. We will ini- 
tiate here the latter, by understanding how small vac- 
uum fluctuations behave in bimetric theories and mas- 
sive gravity. As a by-product, we are able to understand 
stability properties of BHs in these theories and begin 
to understand how gravitational waveforms differ from 
general relativity [see also Ref. [58( for a recent attempt] . 



Massive gravitons and the Gregory-Laflamme 
instability 



1 Note however that the theory considered in Ref [54j does not 
satisfy the Fierz-Pauli tuning and hence it contains a ghost. It 
would be interesting to repeat such calculation for viable theories. 
In this case however, the Vainshtein mechanism discussed in the 
main text may prevent a consistent linear analysis. 



In a very recent paper [241) I. Babichev and Fabbri 
showed that the mass term for the graviton can be 
interpreted as a Kaluza-Klein momentum of a four- 
dimensional Schwarzschild BH extended into a fiat higher 
dimensional spacctime. Such "black string" spacctimes 
are known to be unstable against long-wavelength per- 
turbations, or in other words, against low-mass pertur- 
bations, which are spherically symmetric on the four- 
dimensional subspace. This is known as the Gregory- 
Laflamme instability [5 3 . |60| , which in turn is the analog 
of a Raylcigh-Plateau instability of fluids [61], [62| . Based 
on these results, Ref. J2J] pointed out that massive tensor 
perturbations on a Schwarzschild BH in massive grav- 
ity and bimetric theories would generically give rise to a 
(spherically symmetric) instability. In the following we 
confirm these results within a more generic framework 
and extend them to generic modes and to the case of 
Schwarzschild-de Sitter BHs. 

One of the important open questions is the end-state of 
such instability. For black strings, there is reasonable ev- 
idence that break-up occurs [63(. But the spacetimes we 
deal with are spherically symmetric, and so is the unsta- 
ble mode. A possible end-state is a spherically symmetric 
BH endowed with a graviton cloud (see e.g. Ref. J64|). 
An analysis of the nonlinear equations in case of spherical 
symmetry is left for future work. 



D. Massive bosons and BH superradiance 

The interaction of generic bosonic fields with spinning 
BHs gives rise to interesting phenomena, related to BH 
superradiance [18M21J . Due to the dissipative nature of 
the BH horizon and to the existence of negative-energy 
states in the ergoregion of a spinning BH, low-frequency 
w monochromatic bosonic waves scattered off rotating 
BHs arc amplified whenever the following condition is 
met, 



uj < mfl 



H 



(1) 



where f2# is the angular velocity of the BH horizon and m 
is an integer characterizing the azimuthal dependence of 
the wave. The extra energy deposited in the wavepacket's 
amplitude is extracted from the BH, which spins down. 

Superradiance is prone to very interesting "side- 
effects," such as BH bombs [20. 
[m,[6i,[62i and BH instabilities [ll|, 
review see Ref. [75(). 

The amount of energy extracted through superradi- 
ance strongly depends on the spin of the field. Mass- 
less spin-2 (gravitational) waves can be amplified ~ 300 



.floating orbits 
lH,|6i-|7| (for a 



Ref. 1241 appeared while our work was on its last stages. 



times more than scalar waves. Superradiance scattering 
for massive waves with nonvanishing spin is much more 
involved, due to spin-spin coupling effects [n|. How- 
ever, a generic expectation is that superradiant instabil- 
ities triggered by massive bosons are more effective for 
higher spin. Finally, even in the scalar case superradiant 
effects might be enormously amplified due to the inter- 
action with ordinary matter |76[. 

We are particularly interested in supcrradiance- 
triggcrcd BH instabilities which are sustained by mas- 
sive fields. Ultralight bosons have received widespread 
attention recently as they are found in several exten- 
sions of the Standard Model, for instance in the string 
axiverse scenario [f|, Q where a plethora of massive 
pseudo-scalar fields called axions covers each decade of 
mass range down to the Hubble scale and fields with 
10 _22 eV < m s < I0 _10 eV arc of particular interest 
for BH physics 77(j. In parallel, massive hidden U(l) 
vector fields also arise in extensions of the standard 
model @, y, Hj, UM i highlighting the importance of un- 
derstanding the physics of such fields around BHs. 

Superradiant instabilities were studied extensively for 
scalar fields both in the fre quen cy- and in the time- 
domain 0, [M [H EH [HI [ZSIZI] ■ The non-separability 
of the field equations for a massive vector field in a Kerr 
background has hampered its study for decades (see for 
instance Ref. [79| for some references on the nonrotating 
case). Very recently however, progress has been made. In 
the frequency domain slow-rotating expansions were used 
to prove that massive vectors are superradiantly unsta- 
ble [f0|, OH , these results were confirmed using evolutions 
of wavepackets around Kerr BHs [30[ . It was shown that 
the massive vector field instability can be orders of mag- 
nitude stronger than the massive scalar field. 

The instability is regulated by two parameters, the BH 
spin a/M and the dimensionless parameter M\i (in units 
G = c = 1), where M is the BH mass and m g = fih 
is the bosonic field mass. For ultralight scalar fields 
around massive BHs, the instability timescale can be 
of the order of seconds for solar-mass BHs and of the 
order of hundreds years for a supcrmassive BH with 
A/ ~ 10 9 M Q [l|, [s|, [33|j typically much shorter than the 
evolution timescale of astrophysical objects. The insta- 
bility timescale for spin-1 massive fields can be up to 
three orders of magnitude shorter [T3,[If],[3!|. To summa- 
rize, this mechanism can be very efficient for extraction of 
angular momentum away from the BH. As a consequence, 
observations of massive spinning BHs can effectively be 
used to impose bounds on ultralight boson masses |10| . 



E. Framework 

We wish to describe two different cases: i) the interac- 
tion of a generic massive spin-2 field with standard grav- 
ity, that is, we consider the massive tensor as a probe 
field propagating on a geometry which solves Einstein 
equations; ii) the linearized dynamics of a massive gravi- 



ton as it emerges in nonlinear massive gravity. It turns 
out that both cases can be described consistently within 
a common framework. 

More specifically, we consider the action for two tensor 
fields, g^ and / M „, with a ghost- free nonlinear interac- 
tion between them (cf. Eq. ([5]) below). This class of 
theories is usually referred to as "bimetric gravity" [fa - 
\u\ . The fluctuations of the two dynamical metrics can 
be separated and describe two interacting gravitons, one 
massive and one massless. 

Nonlinear massive gravity [fj, uM, LL21 is obtained from 
the bimetric theory in the limit where the field f^ v be- 
comes nondynamical, i.e. taking M/ — >• in Eq . © 
and considering /^„ as a given auxiliary field [50|. In 
this limit, / M „ can be interpreted as a background metric 

in which the linearized massive fluctuations nf^ propa- 
gates. On the other hand, g^ is a solution of the full non- 



,(m) 



linear field equations such that we have g^ v = /,, 

A crucial point is to identify the background solution 
over which the massive tensor perturbations propagate. 
Linearization of massive gravity is typically considered 
around a flat, Minkowski background. Here instead we 
wish to describe the linearized dynamics around a non- 
linear vacuum solution, i.e. a BH geometry. Regular, 
nonlinear, solutions in bimetric and massive gravity are 
challenging to find and they might exhibit a rich struc- 
ture [80N83J . In bimetric theories new curvature invari- 
ants, such as / = g^ffiv, can become singular at the 
horizon. It was shown that the only way to avoid a 
singular horizon is to req uire both metrics to have co- 
incident horizons |80l . |82|. The same arguments were 
used to show that regular BHs can exist in massive grav- 
ity theories with a flat nondynamical metric provided at 
least one of the metrics is non-diagonal (or non-stationary 
and axisymmetric) when written in the same coordinate 
patch [801 ] . 

In massive gravity the diffeomorphism of general rel- 
ativity is broken, so in principle one is not allowed to 
change coordinates to avoid this problem. This im- 
plies that, assuming a flat background, BH solutions in 
Schwarzschild coordinates must have a component g tT to 
avoid a singular horizon. This component implies a time- 
dependence and nonzero energy flux T tr near the horizon, 
which might even lead to the disappearance of BHs in this 
theory [84J]. Due to the Yukawa-like potential the BH 
gravitational field is screened by a negative energy den- 
sity which is accreted by the BH because of the ingoing 
flux T tr leading to a decrease of the BH mass. Although 
the timescale should be much longer than the Hubble 
time (and hence astrophysically irrelevant), it seems to 
be an anomaly of massive gravity. 

To avoid dealing with such problems, we consider the 
special case in which the background solutions are the 
same as in general relativity. In bimetric theories this 
can be accomplished by taking the two metric to be pro- 
portional, fuu = C 2 g^ v , as discussed in detail below [see 
also Ref. [49|]. This choice also avoids the singular hori- 
zon problem, as the two metrics have the same horizon. 



The linearized equations describing the fluctuations of 
the two metrics can be easily decoupled and they de- 
scribe one massless graviton (which is described by usual 
linearized Einstein dynamics), and a massive graviton 
which is described by the Fierz-Pauli theory on a curved 
background [HlsEj]. 

On the other hand, in the limit of massive gravity this 
is equivalent of taking the nondynamical metric as being 
the BH spacetime instead of the usual flat spacctime. Al- 
though perfectly consistent with the field equations, this 
choice seems somewhat unnatural and other nonlinear 
background metrics can be considered [cf. Ref. [85[ for 
a recent review] . The fluctuations of the physical metric 
Qiiv propagate on a nonlinear BH background / M „ and 
they are also described by Fierz-Pauli theory. 



For massive spin- 2 particles wc must have 2s + 1 = 5 
degrees of freedom. The only choice for the constant k 
that describes a single massive graviton is the Fierz-Pauli 
tuning, k = 1 [22[. In this case, the full set of linearized 
equations reads: 



(□- A t 2 )V=0, 0"V = °> 



0. 



(7) 



On the other hand, for k ^ 1 the theory propagates 6 
degrees of freedom. The extra polarization comes from 
a scalar ghost (a scalar with negative kinetic energy) of 



mass m 



— 2(i-k) ^ 2 ^ wmc h arises from the trace 



equation ([5]). The ghost mass approaches infinity as the 
Fierz-Pauli tuning is approached, so that the ghost de- 
couples in this limit. 



III. LINEARIZED MASSIVE GRAVITY ON 
CURVED SPACETIME 

A. The Fierz-Pauli tuning in flat spacetime 

Let us start by reviewing the classical Fierz-Pauli the- 
ory describing a massive spin-2 field in four-dimensional 
flat spacetime. The action is given by [221 



1 



■>FP 



d 4 x 



1 



d x h^d x h^ + d,h vX d v h^ 



16ttG„ 
-d^Wdvh + ^-d x hd x h - tL. (h^h» u - nh 2 ) 



where h = rf^h^ is the trace of the symmetric tensor 
field h^v, if 11 ' is the Minkowski metric, k is an arbitrary 
constant, and fi is the graviton mass. When /i = 0, the 
action reduces to the linearized Einstein-Hilbert action. 
When fi j£ 0, the mass term violates the diffcomorphism 
invariance of general relativity, i.e., this action is not in- 
variant under infinitesimal transformations of the form 



Sh^ = d^ v {x) + d^ix) 



(2) 



The equations of motion are given by [see Ref. |23j for 

a review! 



5S 



DK 



d\d»h* - d\d v h x + r\^ v d\d a h 



\o 



+ dfj,d y h - rj^Dh - /i 2 (h^ - nrj^h) = . (3) 

Acting with d^ on ([3]) we find the constraint 

d v h vti - nd^h = . (4) 

Note that for k = 1/2 this corresponds to the harmonic 
gauge in linearized general relativity Plugging this back 
into the field equations and taking the trace, wc find 



2(1 - n)Dh + (1 - 4k)^ 2 /i = . 



(5) 



Substituting the trace condition, Eqs. (J3J) reads 

(□ - /i 2 ) V = (2k - 1) [d„d v h + v^^h] ■ (6) 



B. Massive spin-2 particles on curved spacetimes 

Let us now generalize the equations of motion for mas- 
sive spin-2 particles on a curved background |49l 18a. l87j] . 
The more general ghost-free action of two interacting 
spin-2 fields, without matter couplings, is given by [17| 



5* = 




J -R f -2M$V{g,f) 



(8) 
where R g and Rf are the Ricci scalars corresponding to 

g^ v and / M „, respectively; M g 2 = 167rG, MJ 2 = 16-kQ 
are the corresponding gravitational couplings, and M v is 
written in terms of M g , Mf and of the parameters of the 
potential term. The quantities /, g denote the determi- 
nant of the respective metric. There is a unique prescrip- 
tion for the latter in terms of only five interaction terms 
which is free from the Boulwarc-Dcscr ghosts on generic 
backgrounds. We schematically denote the potential as 



V 



4 

E«W. -f V =(y/F T f) u (9) 



n=0 



where /3j are coupling constants. The precise form of 
the potentials V„, is not crucial here and we refer to the 
original papers [lj, [Tj| [l7| • 

Although the action © describes a vacuum bimetric 
theory, it reduces to massive gravity in the limit Mf —> 0, 
in which case the kinetic term of the metric f^ vanishes 
and the field is taken to be auxiliary |50| . 

From the action (|5|) we find two sets of Einsteins equa- 
tions for g^ and / M „ 

RMg)-\9^R(g) + ^T^(j) = o, (10) 

1 A/f4 

R,AI) ^UuR(f) + jfiVAi) = o , (ii) 

where the "graviton" stress-energy tensors T.?„ and Tl v 
depend on 7 M v and are defined, e.g., in Ref. [49| . 



Since we want to consider a BH geometry as back- 
ground, we first need to find a BH solution of the field 
equations. As previously discussed, this is a challenging 
and controversial issue [see Ref. [85| for a recent survey 
of hairy BHs in massive gravity] . 

Here we make the simplest choice and consider two 
proportional background metrics f^ = C g^ v (we use 
the bar notation to denote background quantities). Re- 
markably, in this case the solutions coincide with those 
of general relativity. Indeed, Eqs. (|10|) and ([TTj) reduce 

to Hi 



Rpv — 7}9pvR + A g <y M „ — , 



(12) 



which are just two copies of Einstein's equations with 
two different cosmological constants. The latter are writ- 
ten in terms of the parameters of the interaction poten- 
tials and of the gravitational couplings [49|. Further- 
more, consistency of the background equations requires 
h. g = Af, which translates into a quartic algebraic equa- 
tion for the constant C. Classical no-hair theorems of 
general relativity guarantee that the most general sta- 
tionary BH solution in vacuum and with a cosmological 
constant is the Kerr-(Anti) de Sitter metric. Therefore, 
when A g = Af > the fields g^ u and / M „ describe two 
identical Kerr-dc Sitter BHs. 

Since we are interested in local physics near massive 
BHs, we shall consider A g w ps Af. This condition 
can be satisfied exactly by requiring a fine tuning of 
the interaction couplings [49(. Alternatively, even with- 
out fine tuning, realistic values of the cosmological con- 
stant should not play any role in describing local physics 
at the scale of astrophysical compact objects. There- 
fore, we can safely neglect those terms and focus on 
asymptotically-flat Kerr BHs as background solutions. 
In Boyer-Lindquist coordinates, these are described by 
the line element: 
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Note that the perturbations are generically independent, 
$Qnv 7^ fi.ffj.u- From Eqs. (fTU)) - ([TTj) . the linearized field 
equations read 
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and £ p 1 is the operator representing the linearized Ein- 
stein equations in curved spacetimes: 
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where we already assumed A g = = Af. 

Taking appropriate linear combinations of the metric 
fluctuations, 
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linear equations decouple: 
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(21) 

(22) 
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From the equations above, it is clear that the theory de- 
scribes two spin-2 fields, hpj and hfX . The former is 
massless and it is described by the linearized Einstein- 
Hilbcrt action, whereas the latter has a Fierz-Pauli mass 
term defined as 
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(24) 



+ M 2 b? cos 2 9, A = (r-r + )(r-r-), r± 



where E 

M(l ± y/1 — a) and a = J/M 2 . This spacetime describes 
a rotating BH with mass M and angular momentum J 
in G = c = 1 units. 

Let us now consider fluctuations around the back- 
ground metrics: 
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Note that not all parameter /3j in the equations above are 
independent [46j |. 

What we have discussed so far is valid for bimetric 
theories (|SJ). It is worth stressing that linearized massive 
gravity can be recovered taking the limit Sf^ — > and 
M f -> in Eq. flU]) such that 5fa/Mj -> 0. In this 
limit only Eq. (|16p survives as dynamical equation. In 
the massive gravity limit, this equation can be written in 
the same form as in Eq. (|23j) for the perturbation bgp V , 
but with a mass term 
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(25) 



Therefore, also in this case the theory describes a massive 
graviton propagating in the curved background g^ v = 

U v /c 2 . 

We have just proved that in both cases (bimetric the- 
ories and massive gravity) the linearized equations de- 
scribing a massive spin-2 field on a curved spacetime are 
described by an equation of the form (|23|) . In the case of 
bimetric theory one also has Eq. (|22[) , which we ignore 
since it describes a standard massless graviton and it is 
decoupled. 

In flat spacetime, the equations of motion (f2"3")l reduce 
to Eq. ^ whereas, on curved background they reduce to 
the system: 
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(27) 
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where, here and in the following, we have suppressed the 
superscript "(m)" for simplicity. This set of equations 
can be shown to be the only one that consistently de- 
scribes a massive spin-2 coupled to gravity in generic 
backgrounds |87j . In the rest of this paper we will in- 
vestigate Eqs. (|2l))) - (f2"g|) on a BH background. 



IV. INSTABILITY OF BLACK HOLES 

AGAINST SPHERICALLY SYMMETRIC 

FLUCTUATIONS 

We start by showing that Schwarzschild BHs are gener- 
ically unstable against spherically symmetric perturba- 
tions |24| . This is a generic and strong instability, as 
we will show. To lay the necessary framework, con- 
sider a generic tensor field h^ in a Schwarzschild back- 
ground. Due to spherical symmetry, the tensor field h^ u 
can be conveniently decomposed in a complete basis of 
tensor spherical harmonics [89|, |9CJ. Furthermore, the 
perturbation variables are classified as "polar" or "ax- 
ial" depending on how they transform under parity in- 
version (9 — > it — 6, 4> —> 4> + 7r )- Polar perturbations 
are multiplied by ( — 1)' whereas axial perturbations pick 
up the opposite sign (— 1)' +1 . We refer the reader to 
Refs. [25], [9l[ for further terminology used in the litera- 
ture. 

We decompose the spin-2 perturbation in Fourier space 
as follows: 

/+oo 
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where h^ lm and hP°} ar ' lm 
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(29) 



are explicitly given in Ap- 
pendix [A] In a spherically symmetric background, the 
field equations do not depend on the azimuthal number 
m and they are also decoupled for each harmonic index 
I. In addition, perturbations with different harmonic op- 
posite parity decouple from each other. 



The details of the perturbation equations arc provided 
in Appendix [X] In this section, we are only interested in 
the I = polar sector. The perturbations G, r/o and r\\ 
as given in Eq. (|A2[) are not defined for I = because 
their angular dependence is vanishing. The remaining 
dynamical variables can be recast into a simple monopolc 
equation. First, wc use the constraints (|A20[) and (|A17[) 
to eliminate Hq and H2 as defined in Eq. (JA2I) . Then, 
we use a generalization of the Bcrndtson-Zcrilli transfor- 
mations: 
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After substituting these transformations into the system 
of equations we arrive at a single wave equation of the 
form: 



dr 2 



ip a + [u 2 -V (r)] ip o = 0, 



(30) 



with 
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In this form it is clear that in the massless limit the 
mono polc reduces to the scalar-field wave equation with 
I = [1|. 

We have solved Eq. (|3"0j) subjected to appropriate 
boundary conditions (regularity at the horizon and at in- 
finity, see also next sections) by direct integration, look- 
ing for eigenvalues uj = ujr + iuij. Given the time depen- 
dence (|2"9"|) , stable modes are characterized by u>i < and 
unstable modes by u>j > 0. Wc found one unstable mode, 
detailed in Fig. [1] and characterized by a purely imagi- 
nary, positive component. This is a low- mass instability 
which disappears for M[i > 0.43 and has a minimum 
growth timescale of around Mu>i ~ 0.046. In fact, as 
recognized very recently [24| while our own work was in 
its final stages, the linearized equations (|26p are equiv- 
alent to those describing four-dimensional perturbations 
of a five-dimensional black string after a Kaluza-Klein 
reduction of the extra dimension. Therefore, the system 
is affected by Gregory-Laflamme instability j59l . |6C| that 
manifests itself in the spherically symmetric, monopolc 
mode. One interesting aspect of our own formulation is 
that we are able to reduce this instability to the study of 
a very simple wave equation, described by (|30[) . 

To summarize, in this setup Schwarzschild BHs are un- 
stable. The instability timescale depends strongly on the 
mass scale [i. For low masses, we find numerically that 
luj ~ 0.7/i, in good agreement with analytic calculation 
by Camps and Emparan [62|. 

The Gregory-Laflamme instability only affects 
spherically-symmetric (I = 0) modes [601 ]. so we expect 
the rest of the sector to be stable. We confirm this result 
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FIG. 1. Details of the instability of Schwarzschild (de Sitter) BHs against spherically symmetric polar modes of a massive 
spin-2 field. The left panel shows the inverse of the instability timescale uji = 1/r as a function of the graviton mass /i for 
different values of the cosmological constant A 9 = A/, including the asymptotically flat case A 9 = 0. Curves are truncated 
when the Higuchi bound is reached fj 2 = 2A 9 /3 [88|| . For any value of A 9 , unstable modes exist in the range < M\i < 0.47, 
the upper bound being only mildly sensitive to A 9 . The right panel shows some eigenfunctions in the asymptotically flat case. 
The eigenfunctions decay exponentially at spatial infinity and are progressively peaked closer and closer to the BH horizon for 
masses close to the threshold mass M/x ~ 0.43. 



in Sec. [V] below, where we derive the complete linear 
dynamics on a Schwarzschild metric. 

A more relevant question is related to the role of a cos- 
mological constant. When the background metrics are 
two copies of Schwarzschild-de Sitter solutions, the field 
equations (|26p do not arise from a Kaluza-Klein decom- 
position of a five-dimensional black string. Thus, it is 
not obvious a priori if the monopole instability discussed 
above survives when A g = Kf ^ 0. 

Our formalism can be immediately extended to ac- 
commodate Schwarzschild-de Sitter backgrounds. In this 
case, Eq. (f2U)l is modified with new terms proportional 
to A g , see e.g. Eq. (2.1) in Ref. [92(. From the latter 
equation, one obtain the same divergenceless and trace- 
less conditions as in Eqs. (|27|) and (|28[) . Finally, using 
these conditions and the commutator of two covariant 
derivatives, it turns out that the linearized field equation 
is precisely as in Eq. (|26p . That is, terms that explicitly 
depend on A g cancel out and the only contribution of the 
cosmological constant arises through background quan- 
tities. From the system (|2l))) - (f2"5)) . it is straightforward 
to obtain a master equation for spherical perturbations 
of Schwarzschild-de Sitter BHs. Here we omit the details 
and only give the final result. The monopole is described 
by an equation of the same form as Eq. p0[) . but where 
the potential now reads: 
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Using the same technique as before, we have integrated 
Eq. (f5U|) with the potential (|3~Tj) . The results arc shown 



in Fig. Q]for various values of A g = Af. Note that mas- 
sive spin-2 perturbations propagating in an asymptot- 
ically de Sitter spacetime are subjected to the bound 
fi 2 > 2A g /3 [8f|. Below such bound, the helicity-0 com- 
ponent of the massive graviton becomes a ghost. When 
the bound is saturated, fj, 2 = 2A g /3, the helicity-0 mode 
becomes pure gauge and the instability disappears. The- 
ories with such fine-tuning are called "partially masslcss 
gravities" [§1 [H [see also Refs. [H, m, [H]] and they 
arc not affected by the monopole instability discussed 
above. Finally, as shown in Fig. [1] the instability is even 
more effective for Schwarzschild-de Sitter BHs and it ex- 
ists roughly in the same range of graviton mass. 

For both Schwarzschild and Schwarzschild-de Sitter 
BHs, the instability timescale is of the order of the Hub- 
ble time when (i~ 2x 10~ 33 eV [Hj]. This of course, 
does not mean that the observation of compact objects 
imposes constraints on the graviton mass p. Rather, it 
suggests that the background solution used to describe 
these geometries is likely not the physical one. It would 
seem that a suitable background geometry is given by the 
end-state of this monopole instability. 

Our linear analysis cannot handle the nonlinear devel- 
opment of the instability, nor the nonlinear final state. 
However, from the mode profile in Fig. [TJ it is tempting 
to conjecture that a Schwarzschild BH surrounded by a 
graviton cloud could be a possible solution of the field 
equations. We note that this endstate is completely dif- 
ferent, as it must, from the standard Gregory-Laflamme 
instability which acts to fragment black strings j6l|, [63[ . 



3 Unlike the observation of rotating compact BHs, discussed later 
on, which does impose strict limits on the graviton mass. 



V. MASSIVE SPIN-2 FIELDS ON A 
SCHWARZSCHILD BACKGROUND 

We have established the instability of spherically sym- 
metric fluctuations in non-rotating backgrounds. We now 
generalize the analysis to the full set of non-axisymmetric 
polar and axial perturbations. 



A. Axial sector 

The axial field equations are derived in Appendix [S] 
The axial sector is fully described by the following sys- 
tem: 



equation 
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where A = 1(1 + 1) and we have defined the tortoise co- 
ordinate r* via dr / 'dr„ = / = 1 — 2AI/r. The functions 
Q(r) = f(r)hi and Z(r) = hs/r are combinations of the 
axial perturbations as defined in Eq. (|A1[) . whereas the 
source terms are given by 






Axial dipole mode 



(34) 
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The / = monopolc mode does not exist in the axial 
sector since the angular part of the axial perturbations 
(|Alj) vanishes for I = 0. For the dipole mode (1=1 or 
equivalently A = 2), the angular functions Wi m and Xi m 
vanish and one is left with a single decoupled equation: 
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where s = 0,1,2 for scalar, vectorial, or tensorial per- 
turbations. These transformations were first found by 



Bcrndtson [97J when studying the massless graviton per- 
turbations of the Schwarzschild metric in the harmonic 
gauge. In the massless limit the vectorial degree of free- 
dom can be removed by a gauge transformation, but for 
/i 7^ it becomes a physical mode. Note that the wave 
equation (|37j) for s = 1 is identical to that describing 
electromagnetic perturbations of Schwarzschild BHs |25| > 
thus the axial spectrum of massive spin-2 perturbations 
should include a mode which approaches that of an elec- 
tromagnetic mode in the low- mass limit. 



B. Polar sector 

The polar equations are more involved and derived in 
Appendix (fA). The polar sector is fully described by a 
system of three coupled ordinary differential equations: 
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where the dynamical variables AT, m an d G are defined 
in Eq. (|A2[) and the source terms are given by 

S K = A7^ + 8iKm + A(A - 2)<r 1 ^ + A(A - 2)p 1 G , 

dr 
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2. Axial massless limit 

It is interesting to note that in the massless limit we 
can use the transformations 
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to reduce the system to a pair of decoupled equations, 
given by a "vectorial" and a "tensorial" Rcggc- Wheeler 



The coefficients dj, 0i, 7,, Si, CTj, pi are radial functions 
which also depend on u and I. These equations arc rather 
lengthy and since their explicit form is not fundamental 
here, we made them available online in Mathematiga 
notebooks I32I1. 



1. Polar dipole mode 

The polar monopole was already investigated in Sec- 
tion IIVI and shown to lead to Gregory-Laflamme-like in- 
stabilities [2J|. We now study the dipole mode. In the 
dipole case, I = 1, A = 2, the radial function G iden- 
tically vanishes and we are left with a pair of coupled 
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equations satisfying the following system 



1. Quasinormal modes 



2 d 2 K „ dK - d?7i f 



r?d 2 rii A dni ~ , d/v ? r ^ 



2. Polar massless limit 

In the massless limit we can use the argument pre- 
sented by Bcrndtson in Ref. [97| to reduce the system 
to three decoupled equations, one "scalar", one "vecto- 
rial" (|37|) and one "tensorial" equation described by Zer- 
illi's equation [98[ Q In the massless limit the scalar and 
the vectorial degrees of freedom can be removed by a 
gauge transformation but, for /x ^ 0, they become phys- 
ical. Thus, we expect that the small-mass limit of mas- 
sive gravity spectrum includes a family of modes which 
are identical to that of a scalar and an electromagnetic 
modes (these modes are discussed in Ref. [25[ and avail- 
able online at |32|). 



C. Results 

We have solved the previous systems of equations sub- 
jected to appropriate boundary conditions, which de- 
fines an eigenvalue problem for the complex frequency 
w = ujr + iwi; this pro blem can be solved using several 
different techniques [25l . |26[ which we detail in Appendix 

m 

In general, the asymptotic behavior of the solution at 
infinity is given by 



*i(r) 



B je 



-ik 



+ Cje 
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where koo = vf 2 — w 2 and, without loss of general- 
ity, we assume Rc(fc oc ) > 0. The spectrum of massive 
perturbations admits two different families of physically- 
motivated modes, which are distinguished according to 
how they behave at spatial infinity. The first family in- 
cludes the standard QNMs, which corresponds to purely 
outgoing waves at infinity, i.e., they are defined by Bj = 
[251 . The second family includes quasibound states, 
defined by Cj = 0. The latter correspond to modes spa- 
tially localized within the vicinity of the BH and that 
decay exponentially at spatial infinity [lfj, 123 ILll UM ■ 



4 Note that in these transformations there are four functions. One 
tensorial, one vectorial, and two scalars. However one the scalar 
functions is simply the trace of h M „ , which vanishes in our case (in 
their notation is the scalar function ipo, not to be confused with 
the scalar function used here). We stress again the importance 
of having a vanishing trace in order to have a correct number of 
degrees of freedom. 
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FIG. 2. QNM frequencies for axial I = 1, 2 modes, for a range 
of field masses M\i — 0, 0.04 . . . , 0.52. Points with largest |w/| 
correspond to /i — > 0. The fundamental mode (n = 0, circles) 
and the first overtones (n = 1, triangles) are shown. In the 
massless limit the "vector" modes have the same QNM fre- 
quency as the electromagnetic field, and the "tensor" modes 
have the same QNM frequency as the massless gravity per- 
turbations. 



The axial QNM frequencies for different values of the 
spin-2 mass are shown in Figured As expected, for / > 2 
one can sensibly group the modes in two families for any 
given I and n. They can be distinguished by their behav- 
ior in the massless limit, the spectrum of the "vector" 
modes reduces to the spectrum of the photon, while the 
"tensor" modes, which are the only physical modes in the 
massless limit, approaches the spectrum of the massless 
gravity perturbations. For the lowest overtones, as the 
mass increases the decay rate decreases to zero, reaching 
a limit where the QNM disappears. This is linked with 
the decreasing height of the effective potential barrier as 
was previously discussed in Ref. [99( . The limiting behav- 
ior, when the damping rate reaches zero are the so-called 
quasiresonant modes, whi ch w ere already shown to oc- 
cur for massive scalar [99|, Il00l | and massive vector [l0l| 
fields. 

Polar QNMs are more challenging to compute, because 
the perturbations equations are lengthy and translate 
into higher-term recurrence relations in a matrix- valued 
continued- fraction method [26[. On the other hand, due 
to the well-known divergent nature of the QNM eigen- 
functions [25|], a direct integration is not well-suited to 
compute these modes precisely. Instead of computing 
these modes, in the following we shall rather focus on 
quasibound states - both in the axial and polar sector - 
which are easier to compute [cf. Appendix [B] and more 
relevant for our discussion. 
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2. 



Quasibound states 



Besides the QNM spectrum, massive fields can also 
be localized in the vicinity of the BH, showing a rich 
spectrum of quasibound states with complex frequencies. 
Here the terminology 'quasi' stands for the fact that these 
states decay due to the absorption by the BH, hence 
the complex frequencies. Bound states were alr eady 
considered for massive scalar [7l|, Dirac J1021 Il03| and 
Proca [zl, Il04( fields. In the small-mass limit Mix < I, 
it was shown that for these fields the spectrum resembles 
that of the hydrogen atom: 



Wfi/jU ~ 1 - 



{Myf 



2(j + 1 + nf 



(46) 



where j = I + S is the total angular momentum of the 
state with spin projections 5 = — s, — s + 1, . . . , s — 1, s. 
Here s is the spin of the field. For a given I and 
n, the total angular momentum j satisfies the quan- 
tum mechanical rules for addition of angular momenta, 
\l~s\<j<l + s. 

Our results show that the spectrum (J46J also describes 
massive spin-2 perturbations which is also confirmed an- 
alytically for the axial mode I = 1 (see Eq. ()D10I) of 
Appendix [D]). In Fig. [3] we show the quasibound-state 
frequency spectrum for the lowest modes. Apart from 
the polar dipole (we discuss this in detail below), all 
other modes follow a hydrogenic spectrum as predicted 
by Eq. (|4^|) . The monopole I = [which belongs to a 
different family than the unstable monopole mode dis- 
cussed in Sec. |IV] is fully consistent with 5 = +2 which 
is in agreement with the rules for the sum of angular mo- 
menta, \l — s\ < j < I + s =>■ j = 2. For each pair 
I > 2 and n there are five kinds of modes, characterized 
by their spin projections. Here we do not show the mode 
I = 2, n = 0, S = 1, which is very difficult to find numer- 
ically due the complicated form of the polar equations 
and his tiny imaginary part. Besides that, the existence 
of the mode 1 = 2, n=l, S = with approximately the 
same real frequency makes it even more challenging to 
evaluate the I = 2, n = 0, S = 1 mode with sufficient 
precision. 

Evaluating the dependence of WjQu) in the small-M/i 
limit turns out to be extremely challenging, due to the 
fact that uji is extremely small in this regime. Our re- 
sults indicate a power-law dependence of the kind found 
previously for other massive fields [70|, u>i/fx oc — (M/x) 7 ', 
with 



7] = 4/ + 25 + 5 . 



(47) 



The fact that the modes I = L, S = Si and I = L + S\, 
S = — Si have the same exponent is a further confir- 
mation of this scaling. Note that only the constant of 
proportionality depends on the overtone number n and 
it also generically depends on I and 5. This is confirmed 
analytically for the axial mode I = 1, S = 1 ,n = 0, as 
shown in Fig. 21 where we see that in the low-mass limit 



the numerical results approaches the analytical formula 
derived in Appendix [D] given by 
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(48) 



The quasibound state found for the polar dipole is 
clearly the more interesting. This mode appears to be 
isolated from the rest of the modes and it docs not follow 
the small- mass behavior predicted by Eqs. (|46p and (|47p . 
Furthermore, we have found only a single fundamental 
mode for this state, and no overtones. For this mode, 
the real part is much smaller than the mass of the spin-2 
field. 

The real part of this special mode in region M\i < 0.4 
is very well fitted by 



w R /H w 0.72(1 - Mix) 



(49) 



For the imaginary part we find in the limit Mil <C 1, 

w//A*«-(M/i) 3 . (50) 

That this mode is different is not completely unexpected 
since in the massless limit it becomes unphysical. This 
peculiar behavior seems to be the result of a nontrivial 
coupling between the states with spin projection 5 = — 1 
and 5 = 0. Besides that, this mode has the largest bind- 
ing energy {ujr/ li— 1) for all couplings M/j,, much higher 
than the ground states of the scalar, Dirac and vector 
fields (see Fig. 7 of Ref. [79]). However the decay rate is 
very large even for small couplings M\x, corresponding to 
a very short lifetime for this state. 

To summarize, the I > modes of Schwarzschild BHs 
in massive gravity theories arc stable, with a rich and 
potentially interesting fluctuation spectrum, which could 
give rise to very long-lived clouds of tensor hair in the 
right circumstances. We now show that once rotation 
is included, this hair grows exponentially and extracts 
angular momentum away from the BH. Thus, while the 
monopole / = mode is unstable even in the static case, 
the / > modes suffer for a supcrradiant instability only 
above a certain threshold of the BH angular momentum. 



VI. MASSIVE SPIN-2 PERTURBATIONS OF 
SLOWLY ROTATING KERR BHS 

In Ref. [1J| a method to study generic perturbations of 
slowly rotating BHs was developed. Here we extend this 
method to massive spin-2 perturbations of slowly rotat- 
ing Kerr BHs. We derive the linearized field equations to 
first order in a, although our analysis can be generalized 
to higher order in the BH angular momentum. 

The technique is detailed in Appendix [C] and it con- 
sists in a decomposition of the perturbation equations in 
tensor spherical harmonics and in a expansion in the BH 
angular momentum. The method was originally devel- 
oped to study the gravitational perturbations of slowly- 
rotating stars [27H29| and it has been recently applied 
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FIG. 4. Comparison between the numerical and analytical 
results for the the axial mode 2 = 1, S — 1, n = as a function 
of the mass coupling Mfi. The solid line shows the numerical 
data and the dashed shows the analytical formula (|48[) . 



to BH spacetimes [ill, |3l| . As a result of using a basis 
of spherical harmonics in a nonspherical background, the 
perturbation equations display parity-mixing and cou- 
pling among perturbations with different harmonic in- 
dices. However, as discussed in Ref. UJl, to first order in 



a the eigenvalue spectrum is described by two decoupled 
sets, one for the axial and one for the polar perturbations, 
and all harmonic indices decoupled. In the following we 
discuss the axial and polar sector separately. 



A. Axial equations at first order 

The field equations are derived in Appendix [C] where 
the method to separate the equations is shown. By defin- 
ing: 



h 2 {r) = Z{r)r ( 1 - 



~amM 2 (A + 2) 

Ar 3 w 
amM 2 (A - 2) 

Ar 3 u 



(51) 
(52) 



we obtain that a fully consistent solution at first order is 
such that Z and Q satisfy the following equations: 



d 2 Q 
dr 2 
d 2 Z 
dr 2 



VqQ{t) = S Q Z(r) 
- V z Z(r) = S z Q(r) 



(53) 
(54) 
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with 



2 4amM 2 w 
V Q = w 3 



A + 4 161/ o . , r2 6(4r-9M)(A + 2) 
r r^ Ar°uj 



(55) 



Fz=w 2 



AamM 2 oj 



A - 2 2M 



H g- + /i + amM 



,6(A-2)(r-3M) 



Ar 6 c 



(56) 



-amM' 



Sz = 2/ 



-2)/ 


V-3M 






(6M(4 + A) - r (10 + 3A + 3r 2 w 2 )) " 


. (57) 


Ar 6 w 


"l , (-10 + 3A + 3r V) ~ 
— + amM 2 ± r—= l - J - 

r z Ar°LO 




(58) 



These equations reduce to Eqs. (|32p and (|33|) in the non- 
rotating limit. In the dipole case / = 1, A = 2, the func- 
tion iT vanishes and we are left with a single decoupled 
equation: 



d 2 Q 
dr 2 



V Q Q(r)=Q. 



(59) 



Here the horizon angular velocity Cljj = a/(2Mr + ) was 
expanded to first-order in rotation. When kn < an 
obse rver at infinity will see waves emerging from the 
BH [l05l j. This cor responds to the superradiant condi- 
tion ui < mQ,H [1061 ] . which at first-order in the rotation 
amounts to 



a > 



AMuj r 



(62) 



where ujr is the real part of the mode frequency, u = U)r+ 
iuj. All the polar and axial equations can be brought to 
a form such that the near-horizon solution is given by 
Eq. (f6"0|) . We thus expect that superradiance will also 
occur for massive spin- 2 fields even at first-order in the 
rotation. 

Superradiant scattering leads to instabilities of bosonic 
massive fields [10|, bMj, S, [L0, L73| ■ This instability was 
explicitly shown for scalars and vectors, but generic ar- 
guments indicate that it is present for other integer-spin 
fields. Note that with our convention, unstable modes 
correspond to u>i > 0. These superradiant instabilities 
occur only for waves localized in the vicinity of the BH, 
i.e., quasibound states, so we focus on these states in the 
next sections. 

The continued fraction method can be used to deter- 
mine the quasibound state frequencies of the axial equa- 
tions by imposing an appropriate ansatz which in this 
case is given by 



S>,r) = /W^rV'^flW/W 



(63) 



B. Polar equations at first order 

In line with the non-rotating case, for the polar sector 
we obtain at first order in a three coupled equations for 
K, r)i and G, which generalize Eqs. (|55]). ([3"5]). and ([10]) , 
but in this case the coefficients dj, Pi, 7*, Si, &i, pi are 
also functions of ma. Due to the length of the equations 
we do not show them explicitly here but we made them 
available online in Mathematica notebooks 13211. 



C. Superradiance and quasibound states 

Interesting phenomena, such as BH superradiance, are 
already manifest at first order in the BH angular momen- 
tum. A second order approximation would be necessary 
to consistently describe superradiance (see e.g. Ref. |10J ]) 
but this is beyond the scope of this work. 

As for the Schwarzschild case, at the horizon we must 
impose regular boundary conditions, which corresponds 
to purely ingoing waves 



*iM 



-iknr* 



as r* — > — 00, where 



?H 



u 



m£l 



H 



U) 



UI 



+ 0{a 6 ). 



(60) 



(61) 



where v = —q + ui 2 /q. To compute the quasinormal 
mode frequencies we use q = —\/p 2 — oj 2 and for the 
quasibound state frequencies q = \/ p 2 — oj 2 . Inserting 
Eq. (|63|) into Eq. (|59|) leads to a six-term recurrence rela- 
tion which can be reduced to a three-term rec urrence re - 
lation by successive Gaussian elimination steps 107lll08| . 
For I > 2 we find a six-term matrix-valued recurrence 
relation which can also be brought to a three-term re- 
currence relation using a matrix-valued Gaussian elimi- 
nation. The explicit form of the coefficients is not shown 
here for brevity but it is available online [32J. 

Although the continued-fraction method works very 
well for quasibound states, the multiple matrix inversion 
of almost singular matrices (since some matrices are pro- 
portional to a) makes it very difficult to compute the very 
small imaginary part of the axial quasibound states. We 
therefore use the direct integration method for both the 
polar and axial quasibound states which gives more ac- 
curate results in this case, and use the continued-fraction 
method to check the robustness of our results. 



D. Results 

In the top panels of Fig. [5] we show the absolute value 
of the imaginary part as a function of the rotation param- 
eter for the axial modes I = 1, S = 1 and I — 2, S = — 1. 
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Although a second-order approximation would be needed 
to describe the superradiant regime in a self-consistent 
way [ll| , the first-order approximation predicts very well 
the onset of the instability and should give the correct 
order of magnitude of the instability timescale. For ax- 
ial modes the instability is very weak: even in the most 
favorable cases the instability is almost five orders of 
magnitude weaker than that associated to axial Proca 
modes [ljj, [ljj . This also makes it difficult to track nu- 
merically the axial spin-2 modes with sufficient precision. 
For small masses the real part of the frequency is roughly 
independent on the spin. This is supported by analytical 
results for the axial dipole mode, which can be evaluated 
analytically in the small-mass limit at first order in a [cf. 
Appendix |D] . The analytical formula for the imaginary 
part of the fundamental mode reads 



MW; 



40 
19683 



(a-2r +A i)(M^) n 



(64) 



In Fig. [5] we compare the analytical formula with the 
numerical results for the fundamental overtone and mass 
coupling Mfi = 0.05. Although the imaginary part is 
tiny, the agreement is good in the /i — > limit. Near the 
superradiant regime the agreement is only qualitative, 
as expected since the analytical formula is only valid for 
am/{Mfi) < I. 

The bottom panels of Figure[5]show the imaginary part 
as a function of the BH angular momentum for the polar 
dipole I = 1 and the polar mode I = 2, S = — 2. In this 
case the imaginary part of the mode is larger, and these 
modes are easier to evaluate numerically. The instability 
for the mode / = 2, S = — 2 is roughly two orders of mag- 
nitude weaker than the strongest instability of a Proca 
field [Tq] . Once more the polar dipole mode is the most 
interesting case as it has the largest imaginary part, cor- 
responding to an extremely short instability timescale. 
This agrees with the analysis in the nonrotating case of 
Sec. El where we found that the behavior of this mode is 
different from the rest of the spectrum. 

As shown in the bottom panels of Fig. [5] the polar 
dipole mode displays a peculiar behavior in the super- 
radiant regime, where the power-law dependence is in- 
verted, i.e., the instability is stronger for the lowest mass 
coupling M[i. This suggests that extrapolating the first- 
order results to the superradiant case is probably less 
accurate for this mode. This is confirmed by the behav- 
ior of the real-part of the frequency as a function of the 
spin, as shown in Fig. [7] At first-order the eigenfrcqucn- 
cies can be expanded as 



ujr = u>q + amoji + 0(5 ) 



(65) 



where luq is the eigenfrequency in the nonrotating space- 
time and wi is the first-order correction which is an even 
function of m [1JJ . Hence at first-order we would expect 
that the curves for / = m and I = — m are symmetric 
when reflected around the m = curve. For the polar 
dipole this only happens for very small masses. Note 
also that, contrarily to the rest of the spectrum, the real 



part of the polar dipole mode acquires a nonnegligiblc 
dependence on 5, even in the small fi limit. In fact the 
analytical results for the axial dipole suggest that the 
first-order approximation is only valid for am/ (Mur) < 
I. Since in this case Mujr is much smaller that Mfj,, the 
extrapolation to the superradiant regime is less accurate 
in the polar dipole case. Nonetheless, using the exact 
results in the nonrotating case [cf. Sec. [V] and a linear 
extrapolation of the first-order corrections, we estimate 
the following scaling for the imaginary part of the polar 
dipole mode: 



Mui ~ 7 po iar(arn - 2r + u R )(M^) c 



(66) 



where 7 po iar ~ 0(1) and lor is the zeroth order real fre- 
quency given by Eq. (|49[) . This behavior becomes less 
accurate deep inside the superradiant regime. Although 
such extrapolation is extremely rough, a similar estimate 
has been done in the scalar and in the Proca case and it 
turned out to be very accurate [lQ|. In the scalar case a 
fit similar to Eq. (|66[) agrees with exact results (obtained 
sol ving the Klein-Gordon equation on an exact Kerr met- 
ric [7l|) within a few percents; and, in the Proca case, 
it reproduces the results of exact numerical simulations 
(again in the quasiextremal, 5 ~ 0.99 case) within a fac- 
tor two [301 ] . 

In the case at hand, even if Eq. (|66[) eventually turns 
out to be accurate only at the order-of-magnitude level, 
this would anyway mean that spin-2 fields can trigger the 
strongest superradiant instability among other bosonic 
perturbations. The instability timescale is four orders of 
magnitude shorter than the shortest timescale for Proca 
unstable modes [lOj. A second-order analysis would be 
important to confirm this result, but it will also be very 
challenging. A most promising extension is to perform 
a full numerical analysis (along the lines of Ref. [301 ) m 
the case of massive spin-2 fields around highly spinning 
Kerr BHs. 



VII. DISCUSSION 

The advent of new and powerful methods in BH per- 
turbation theory and Numerical Relativity in the past 
few years allows one to finally tackle traditionally com- 
plex problems. Particularly important to beyond-the- 
Standard-Model physics are scenarios where ultralight 
bosonic degrees of freedom are present; simultaneously, 
massive degrees of freedom turn out to be important out- 
side particle physics, in particular several extensions of 
general relativity encompassing massive mediators have 
been proposed. Thus, the study of massive fluctuations 
around BHs is a timely topic. 

Interesting nonlinear completions of the Fierz-Pauli 
theory have recently been put forward [12 - [l4| . While 
it is at this stage too early to claim a consistent theory 
of massive gravitons (these theories or at least certain 
sectors are either pathological |44l . |48| or phenomenolog- 
ically disfavored [45|). any nonlinear theory describing 
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FIG. 5. Absolute value of the imaginary part of the axial and polar quasibound modes as a function of the BH rotation rate 
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bottom panel: polar dipole mode for I — m = 1. Right bottom panel: polar mode I = m — 2, S — —2. For any mode with 
m > 0, the imaginary part crosses the axis and become unstable when the superradiance condition is met. 
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a massive spin- 2 field - including a massive graviton - 
will eventually reduce to Eqs. (|2l)|) - (|28[) in the linearized 
regime. 

Here we have explored the propagation of massive 
tensors in BH backgrounds as described by Eqs. (|26l) 
and shown that they lead to generic instabilities. 



Schwarzschild and Kerr BHs are both unstable against 
linearized monopole perturbations. These are strong, 
small-mass instabilities whose end-state is unknown. 

Schwarzschild BHs also admit a very rich spectrum 
of long-lived stable states. Once rotation is turned on, 
these long-lived states can grow exponentially and ex- 
tract angular momentum away from the BH. Thus Kerr 



16 



BHs are also unstable against a second mechanism: su- 
perradiance. We showed that the instability is triggered 
when the supcrradiant condition is met, thus providing 
one further and strong piece of solid evidence that super- 
radiant instabilities occur for any bosonic massive field. 
The polar gravitational sector is particularly interesting, 
as it displays the shortest instability timcscale among 
other bosonic fields. Our results arc formally only valid 
in the small BH rotation limit, but previous second-order 
calculations for massive vector fields suggest that a first- 
order analysis provides reasonably accurate results even 
beyond its regime of validity. The most crucial point in 
this regard is the functional dependence of the instability 
timescale for the supposedly more unstable polar dipole 
mode, which we estimate to be: 



Tte 



M{MnY 



7 P oiar(a - 2r + U! R ) 



(67) 



This timescale is four orders of magnitude shorter than 
the corresponding Proca field instability [HJ, [ljj . 

It has been shown that BH supcrradiant instabilities 
together with supermassive BH spin measurements can 
be used to impose stringent constraints on the allowed 
mass range of massive fields [TJJ, LLLI ■ The observation 
of spinning BHs implies that the instability timescale is 
larger than typical competing spin-up effects. For super- 
massive BHs a conservative estimate of these timescalcs 
is given by the Salpcter timescale for accretion at the 
Eddington rate, ts ~ 4.5 x 10 7 years. We find that 
the current best bound comes from Fairall 9 [109J . for 
which the polar instability implies a conservative bound 
/i<5x 10 -23 eV. Unlike bounds for hypothetical mas- 
sive photons, which may interact strongly with matter, 
the previous bound should not be strongly affected by 
the presence of accretion disks around BHs, as the cou- 
pling of gravitons and other spin-2 fields to matter is very 
feeble. 

Our work requires extensions and further analysis (in 
particular, the understanding of the time-development of 
the monopole and superradiant instability requires non- 
linear simulations), and should in fact be looked at as 

I 



the first step in a broader program of understanding 
gravitational-wave emission in massive theories of grav- 
ity. 
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Appendix A: Linearized field equations for a spin-2 
field on a Schwarzschild geometry 



A massive spin-2 field propagates five helicity states 
and one cannot impose the same gauge choices that are 
usually imposed in the massless case. In particular, the 
standard Regge- Wheeler gauge [89| is too restrictive for 
a massive spin-2 field. 

In this paper, we have decomposed the spin-2 field in 
terms of axial and polar perturbations and expanded in 
a complete basis of tensor spherical harmonics. Given 
the expansion (|29p . the axial and polar parts are given 
respectively by 



axial, Im 



(w,r,( 



(0 h l m (u } ,r)csc9d 4 ,Yi m (8,^) 
* h[ m (u J ,r)csced 4> Y lm (e, ( j ) ) 

^ ' ' sin 8 



V 






-/4 m ' 



-h l m (iu,r)sm9deY lm (9,(t ) )\ 
-h[ m (uj,r)smedeY lm (e,(f>) 
h l 2 m (uj,r) sin eWi m {B,<j>) 
h l 2 m (u,r)smeX lm (6,cf>) ) 



(Al) 



t polar ,1m 



(u,r,{ 



ff{r)H l m (u:,r)Y lm H[ m (w , r)Y lm T] l m (tJ,r)d e Y lm 

* f(r)- 1 H% n (u,r)Y hn i 1 l { n (u,r)d e Y lm 

r 2 [K lm (uj,r)Y h 
+G lm (cu,r)W lm 



\ 



^ m ( w ,r)a r im \ 

ri{ m (u,r)d^Yi m 

r 2 G lm (u,r)X lm 

r 2 sin 2 B [K lm (u,r)Yi r , 

-G lm (cu,r)W lm ] 



(A2) 



where f(r) = 1 — 2M/r, asterisks represent symmetric 
components, Y\ m = Yi m (6, </>) are the scalar spherical har- 
monics and 

X lm (6, c/>) = 2d,p [d e Y lm - cot 9Y lm ] , (A3) 

W lm (0, 4>) = d 2 g Y lm - cot 9d e Y lm - esc 2 9d 2 Y lm .(A4) 



1. Axial equations 

The field equations for the axial sector are obtained by 
using the decomposition (|A1|) in Eq. (|26p . Substituting 
into the linearized field equations, we obtain: 



fh 



w 2 -/U 2 + 



A AM 



o 
2-^-/11 = 0, 



/*(, 



/X + -J'" 



r 

2Miw 
'r(r-2M) 



/'i 



2 ,, 2 A + 4 8M 



2(2 -A)/, 
fto+ l , -% =0; 



/X 



2/(r - 3M) 2/ 2 

2 2 "" T 

A - 4 8M 



2 ^i 

r 



W - / ,j 2 + 



h 2 = 0, 



(A5) 
(A6) 
(A7) 
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where A = 1(1 + 1) and / = f(r). Equations (|A5|) and 
(|A6j) correspond to the (tff) and the (r0) component of 
the field equations respectively, and (|A7|) corresponds to 
the (99) component. The transverse constraint (|27p leads 
to the radial equation 



fh[- 2 -^f^h 1 + Jh Q + ^h 2 =0, (A8) 



which can be obtained either from the 9 or the </> compo- 
nent. For the axial terms the trace (|28p vanishes identi- 
cally, 



^axial = _ 



(A9) 



Using the constraint (|A8|) we can reduce the system to 
a pair of coupled differential equations. Eliminating ho, 
we finally obtain the system (|3"2"j) - ([3"3")) . 



2. Polar equations 

Using the decomposition (|A2|) in Eq. (|2l)j) and substi- 
tuting into the linearized field equations, we obtain: 
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(A12) 
(A13) 
(A14) 
(A15) 
(A16) 



Equations (|A10I) - (|A14[) correspond to the equations, respectively. From the (9c/)) component we 
(tt), (tr), (rr), (t6) and (r6) components of the field 
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get Eq. (|A15[) . which combined with the (66) component 
yields Eq. (|A"I5)) . 

The transverse constraint (|27| leads to the following 
radial equations 



fH[ 



2(M -r) 



Hi + iujHq 



A 
— Vo 



0, 



(A17) 



a. Continued-fraction method 

The use of the continued fraction method requires a 
suitable ansatz, which we take to be 

*,-(«, r) - /(r)- 2iMa, r^-^^a«/W n , (B3) 



!H' 2 



2r-3M M 2/ /A 

^ #2 + iioHi + —H K 3-771 = , 

(A18) 

fv'i - 2{M ~ r \ i + ^o + K - (A - 2)G = , (A19) 

for the t, r and component of the constraint, respec- 
tively. Finally, in the polar case the traceless con- 
straint ((28)) yields 



H =H 2 + 2K . 



(A20) 



Unlike the axial sector, the polar equations are not so 
straightforward to further reduce. For I > 2 one could 
use the constraint equations to eliminate 770, ?7i, Hq and 
G and obtain three second-order equations for K, Hi 
and H 2 . However, this choice is not particularly useful, 
because the system do not directly contain the monopole 
and dipole cases (I = 0, 1). For this reason we chose to 
work with K, r\i and G as dynamical variable instead. 

After some tedious algebra, we obtain that the polar 
sector is fully described by Eqs. (|38)) - (j40)) in the main 
text. 



Appendix B: Eigenvalue problem: Quasinormal 
modes and quasibound states 

This appendix details the numerical computation of 
BH eigenfrequencies for massive perturbations. To have 
a well-defined problem we need to define boundary con- 
ditions, and these determine an eigenvalue problem for 
the frequency w, which can be solved using several differ- 
ent tools [25], [26[ . At the horizon we must impose regular 
boundary conditions, which corresponds to purely ingo- 
ing waves 



*iM 



J = 1,2, 



10, 



(Bl) 



as r* — > —00, where $j(r) is any of the radial functions. 
On the other hand, the asymptotic behavior of the solu- 
tion at infinity is given by 



*,-(»•) ~B,-e 



— ik 



M(y.- i ~2u> z ) 



- + C 7 e lfc ~ r : 



MjuT -2u; z ) 



(B2) 

where k^ = y^fi 2 — uj 2 . such that Rc(/c oc ) > 0. For 
massive fields we have to consider two kinds of modes: 
(i) the quasinormal modes (QNM), which corresponds to 
purely outgoing waves at infinity, i.e., they are defined 
by Bj = 0; (ii) quasibound states, defined by Cj = 
and correspond to modes spatially localized within the 
vicinity of the BH and that decay exponentially at spatial 
infinity. 



where v and q are defined as below Eq. (1031) . 



b. Axial dipole 

Inserting (|B3I) into (|3l)|) leads to a three-term recur- 
rence relation of the form 

a ai + f3 a = , 
a„a„+i +/3„a„ +7„a„_i = 0, n>0, (B4) 



where, 

a„ = (n + l)(n + 1 — 4zw) . 



(B5) 



Pn 



-2(n 2 + n-l) 



w 2 (2n-4:iLU + 1) 



- 3q(2n - Aiu + 1) + Ai(2n + 1)uj - Aq 2 + 12w 2 , 



7„ = q 2 (nq + q 2 — 3q — 2iqu) — 
x (nq + q + 3q — 2iqui — u ) 



(B6) 



(B7) 



The QNM or quasibound-state frequencies can be ob- 
tained solving numerically the continued fraction equa- 
tion 



A) 



a 7i 



a ai72 



0. 



(B8) 



This method has been ex tens ively used and described in 
detail elsewhere [ltj, 12a, IllOJ . some routines are freely 
available [32J so we will not discuss it any further. 



c. Axial modes: I > 2 

For I > 2 the axial modes satisfy a pair of coupled 
differential equations, Eqs. (|32[) and (P5|) . Inserting (|B3[) 
into these equations leads to three-term matrix-valued 
recurrence relation, 

aoVi + /3 U = , 
a„U„ +1 + /3„U„ + 7«U„_i = , n > , (B9) 

The quantity U n = I a« , a,\ ) is a two-dimensional vec- 



torial coefficient and ct n , /3 n , 7„ arc 2x2 matrices whose 
form reads, 



a. 



In 



a n 
a r , 



In 6-3A 
7« + 9; ' 



0n = 



Pn A- 2 

-2 p n - 3 / ' 
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with 

0L n = (n + l)(n + 1 — Aiuj) , 

/3„ = 2 - A - 2 (n 2 + n - l) + 



w 2 (2n-4iu; + 1) 



where the superscripts denote a particular vector of the 
(BIO) chosen basis, for example, $J corresponds to bo = 



- 3g(2n - Aiuj + 1) + 4i(2n + l)w - Aq 2 + \2u? , 

(Bll) 

In = q~ 2 [q 2 (n 2 - Ainu - Geo 2 - 9) + 2q 3 (n - 2iu)) 

-2qoj 2 (n~2iLj) + q 4 +uj 4 ] . (B12) 

The matrix-valued three-term recurrence relation can be 
solved using matrix- valued continued fractions [ljj, [Hi ■ 
The QNM or quasibound frequencies are roots of the 
equation MUq = 0, where 



M = A, + a Rl , 



(B13) 



with U n+ i = R^U n and 

K = - (Pn+l + On+iRt+i) 7«+l ■ (B14) 

For nontrivial solutions we then solve numerically 

det|M| = 0. (B15) 

d. Direct integration for quasibound states 

To compute the spectrum of quasibound states a direct 
integration approach is often possible, since the solutions 
asymptotically vanish at spatial infinity, and desirable 
because it converges faster. We start with a series ex- 
pansion close to the horizon of the form 



*i(w,r) 



r^'V^'fr 






rn) 



(B16) 



lO) 



where the coefficients 6„ for n > 1 can be found in terms 
of 6g by solving the near- horizon equations order by or- 
der. We then integrate outward up to infinity where the 
condition Cj = in Eq. (|B2|) is imposed. This allow us to 
obtain the frequency spectrum using a shooting method. 
This method can be extended to solve systems of coupled 
equations [10|, uM ■ Consider a system of N coupled equa- 
tions. Imposing the ingoing wave boundary condition at 
the horizon (|B16[) we may obtain a family of solutions at 
infinity characterized by N parameters, corresponding to 
the iV-dimensional vector of the coefficients bo = {oq }, 
with j = 1, . . . ,N. Note that all the solutions of the 
system of coupled equations must have the form (|B16I) 
near the horizon. We may then compute the bound-state 
spectrum by choosing a suitable orthogonal basis for the 

space of initial coefficients Oq . To do so we perform N 
integrations from the horizon to infinity and construct 
the N x N matrix 



(i) 



(*K 



S m (w) 



lim 



$ 



(2) 

(2) 



v*S> 



v (i) \ 



<o 



(B17) 



{1,0,..., 0}, $j corresponds to b = {0, 1, . . . , 0}, and 

$j ' corresponds to bo = {0, 0, . . . , 1}. The bound-state 
frequency luq = ton + iuj will then correspond to the 
solutions of 



det|5 TO (wo)|=0, 



(B18) 



which in practice corresponds to minimizing det S m in 
the complex plane at arbitrarily large distances. 



Appendix C: Linearized field equations for a spin-2 
field on a slowly rotating Kerr BH 

We will follow Kojima [27| to write the fields equations 
for a spin-2 field in a slowly rotating BH. Since this back- 
ground is still "almost" spherically symmetric we can use 
the decomposition (|29p and insert it in the linearized field 
equations. We can then separate the equations in three 
different groups. 

From the (ti), (tr), (rr), the sum of (89) and (4>4>) com- 
ponents of Eq. (f2l)]) . the t and r components of the trans- 
verse condition (|2T|) . and the traceless condition (|2"8|) . we 
have 



.4 



CO 



m. 



A 2 cos 9 Y lm + B^sm9d e Y 



(/) 



c$d 4l Y lm = o (1 = 0,..., 6) 



(Cl) 



where a sum over (I, m) is implicit, the functions A lr ^ and 
Cj m are some linear combinations of the polar functions 
I/o, Hi,H2, r/o, 771 , K and G. On the other hand ^4 /m and 
_B i?n are some linear combinations of the axial functions 
h , All, h 2 . 

From the (t9), (tip), (r9), (rip) components of Eq. (|2"6]l. 
and the 9, <j> components of Eq. (|2T|) . we have 



a lm + a im COs9)df,Y l 



( C ' (J 

- (/#? +ffi<x*e) (d^/smO) +^(sin0F te ) 
+ ti£* lm + x\2 ( sin ewlm ) = (J = 0,1,2), (C2) 
and 



f/^+Z^cosflW™ 



„( J ) _L X.W 



Hi 
U) y-lr, 



t (J) 



8)(d <p Y lm / S m9)+(^(sm9Y lm ) 



+ x\JX tm - Q„> (sin 9W lm ) = ( J = 0, 1, 2) , (C3) 

where the functions a[ m , /3 Zm , Q m and £ Zm are some lin- 
ear combination of the polar functions, while /3 im , a\ m , 
n im an d Xim belong to the axial sector. 
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From the (9(f>) and the subtraction of 
components of (f2l>|) . we have 

f lm d e Y lm + g lm (d^Y lm /sm6) 
+ (t lm +i lm d^)(W lm / S m6)-- 



and 



0. 



(C4) 



and 



gi m d e Y lm - f lm (d^Y lm / sm6) 
-{t lm + ii m d^){X lm /sm 2 d) 

+ {s lm + s lm d 4> )(W lm /sm9)=Q. 



(C5) 



where fi m , si m and si m are some linear combinations 
of polar functions and gi m , ti m and ti m from the axial 
functions. 

It is easy to see that at zeroth-order in the rotation the 
perturbation equations reduce to 



A {1) - (y {J) 

A lm ~ a lm 



si™ = , (/ = 0, 



,6, J = 0,1,2) 



for the polar sector and to 



4 



(■/) 



th 



0. 



(J = 0,l,2) ; 



(C6) 



(C7) 



for the axial sector, respectively. These equations corre- 
spond to the ones obtained for the Schwarzschild case. 
To separate the angular variables we use the identities 

cos8Y lm = Q l+lm Y l+lm + Q lm Y l ~ lm , 

sm8d e Y lm = Q l+lm lY l+lm - Q lm (l + l)Y [ - lm 



(C8) 



with 



Q 



Im 



(C9) 



(CIO) 



4/2-1 ' 

and the orthogonality properties of scalar, vector and 
tensor harmonics. The separation of the angular depen- 
dence of Einstein's equations for a slowly-rotating star 
was performed in Ref . [27] . Since the above equations are 
formally the same as those considered in Ref. [27|, they 
can be separated in exactly the same way [see Ref. [26| 
for a review]. Below we omit the index m, because in an 
axisymmetric background it is possible to decouple the 
perturbation equations so that all quantities have the 
same value of m. 

From Eq. JUT]) we have [H, H3 



AY'+imCl^+QilA^ + il- 



,(/) 



,co 



i-1 



j(/) 



+ Q l+1 [A$ x -(1 + 2)B l+l 
Equations (JCU) and ([C3} give 



(/ M=0. 



Aa 



(■/) 



;(J) 



Qi(l + l) 



(!-i)(l + 2)t l 

(l-2)(l--t^ (J) 



5(J) 



.(J) 



l)x\-[ +Q- 1)«S 



'ii-i 



U) 



^,( J ) 



S J ) 



Q l+1 1 (l + 2)(l + 3)x\i' 1 -(1 + 2)a^ - ^i 



= 0, 

(C12) 



and 



KB 



(J) 



( J ) , z,( J ) 



(l-l)(l + 2) X r>+a\ J >+r,l 



M) 



t (l + l) (l - 2 )(l -!)$<_[ -(I 



(■/) 



im+4-i 



+ Q l+1 1 (l + 2)(l + 3)^ + (I + 2)/3 ; (J) 



i+i 



AJ) 



= 0. 

(C13) 



Finally, Eqs. flCJ} and JQgJ yield 



A(si +imsi) -imfi - Qi(l + l)gi-i +Qi+\lgi+i = 0, 

(C14) 

A (ti + imii) + imgi - Qi(l + l)/;_i + Qi +1 l f i+1 = . 

(CIS) 
Because the background is nonsphcrically symmetric, 
the radial equations above display mixing between per- 
turbations with opposite parity and different harmonic 
index. To first order, perturbations with given parity 
and harmonic index I arc coupled to perturbations with 
opposite parity and indices J±l. However, as discussed in 
Ref. [ll(, these couplings do not contribute to the eigen- 
value spectrum to first order in a. Finally, neglecting 
the coupling to the opposite parity with harmonic in- 
dices I ± 1, we use Eqs. (|C13|) and (|C15[) to deduce the 
axial equations (|53[) and (|54|) in the main text, while 
the polar equations are obtained from Eqs. ([Clip . (|C12|) 
and (jCT4|) . 



Appendix D: Analytical results for the axial dipole 

In this appendix we generalize Detweiler's analytical 
calculations [33[ for the unstable scalar modes of a Kerr 
BH in the small-mass limit to the case of the massive 
spin-2 axial dipole, to first-order in the rotation. 

Defining R(r) = Q/r the axial dipole equation (|59p 
can be rewritten as 



r 2 f 



dr 



\f 



dR 
dr 



r ui — AamM rui — r 



f(j(j + l) 



2 2 

-(j, r 



2Ms L 



- amM- 



12(4r-9M)\ 



R = 0. 



(Dl) 



where we have defined j = I + S = 2 and s' = 3. From 
now on we consider j and s' to be generic integers and we 
replace their specific values only in the final result (|D14j) 
below. The latter is valid for any j and s' provided j < 
(Gil) s '_ Xo use the method of matching asymptotics we start 
by writing this equation in terms of the dimensionlcss 
variable z = (r — r + )/r + , 



4MV(1 + z) 4 - 2dmMu}(l + z) 



dz \ dz J 
- j(j + 1)Z - 4AfV ^(1 + zf + s' 2 z 
. 3z(l-8z) 



4Mw(l + z) 3 



R = 0. 



(D2) 



21 



where Z = z(z + 1). 

We first expand the equation above forz» 1. For this 
we define the variable x = 4A/fc 00 z and get the equation 

i v ju + iy 

4 



d 2 , s 



.r 
ft, — , 



x 

,2 



a;i? = , 



(D3) 



where we have defined , k 2 ^ = /i 2 — uj 2 , v = M\x - /fcoo and 
have considered u ~ /i. For quasibound states the solu- 
tion of this equation with the correct boundary condition 
at infinity is given by 

- x / 2 x j U(l 



Roo(x) « Cie" 



j-i/,2i + 2,x) 



(D4) 



where Ci is a constant and ?7 (jP, g, x) is one of the conflu- 
ent hypergeometric functions |lll| . For z -C 1, at leading 
order, the behavior of the solution reads 



^oo(r) « C\ 



(Zk^r) 



-(2*o 



v-i-i Hi 



,- r[-i-2j] 
r[-j - 1/] 

2j] 



r[i + j-i/] 



(D5) 



Equation (|D2[) can also be solved in the region where r <C 
max(J/u), j/fJ-)- In this limit, 



4( z fH 



l)Z + s 2 z] R = 0, (D6) 



where we have defined e = 2Af;U, s 2 = s' 2 — ^|-, P = 
—2MJch = -2M(w — mfln) and neglect C(a 2 ) terms in 
P 2 . Note that in order to solve the equation analytically, 

- 2 

we neglect terms C( 22 -), so the approximation is valid 
only if a -C j M/U. 

The solution of the equation above is given in terms 
of hypergeometric functions. Imposing ingoing waves at 
the horizon we get that the general solution is given by 



R H {r) = C 2 e-^(-iy^z^(l + zf 

2-Fi (-3 + iP + a,l+j + iP + a,l + 2iP, -z) , (D7) 

where 2-Fi(a, b, c, z) is the hypergeometric function [11 1| 
and a = \/s 2 — P 2 . Using t he a symptotic properties of 
the hypergeometric function [lllj we can derive the large- 
distance limit z 3> 1 of this solution 



Rh{t) 



C 2 T[1 + 2iP] 

(2M) 1+ J'r[-i 



A)\ 



-j + iP-a]T[-j + iP 
{2M)-*T[1 + 2j] 



(D8) 
T[l+j + iP-a]T[l+j + iP ' 

The near- and far-region solutions have an overlapping 
region when Muj <C j and Mfi <C j and one can find 
a matching condition equating the coefficients of r 3 and 
r —j-i. 

T[2j + l]T[-j - v\ 



r[-2j - i]r[j - v + 1] 
. r[-2j - i]r [j + iP - a 



- 1] T \j + iP + a + 1] 



rpy + i]r [-j + iP - a] r [-j + iP + a] 



(D9) 



At leading order for Mfcoo the right hand side vanishes. 
In the left hand side this corresponds to the poles of 
r [j + 1 — v\ , which are given by iA Q > = j + 1 + n for a non- 
negative integer n, yielding the expected hydrogen-like 
quasibound states. We obtain, to lowest order in M/x, 



MfJL 



j + n+l 



(D10) 



In order to get the imaginary part of the spectrum, we 
expand around this value to get the next-to-leading order 
correction. Writing v = ^ + Sv and assuming 5v <C 1 
we get (for details see e.g. [l!2J ) 






(4fc 0O M) 2 J+ 1 r[-2 J - i]r[2j 



2] 



r[i + 2j] 2 r[2j + 2]r[n + i] 

iP - a + l]T[j + iP + a + 1} 



T[-j + iP-a]T[-j + iP + a} 



(Dll) 



Since there is a pole in one of the T-functions we take to 
lowest order in P and a/e, T[—j + iP — a] w T[—j — s']. 
We then get in this limit 



(-1) 



,_ B , (4fc 0O A/) 2 J+ 1 r[2j 



2]rb' 



2r[l + 2j] 2 r[2j + 2] 2 r[n + l]T[-j + s'} 
x T\j +iP- a + l], (D12) 

where the factor 2 in the denominator comes from a spe- 
cific limit of the T functions and it is related to the fact 
that 1(1+1) is not the exact angular cingenvaluc in a 
rotating background. In the nonrotating limit, the re- 
sult above must be multiplied by a factor 2. [see the 
discussion of Appendix C2 in Ref. [Il[ for details] . The 
imaginary part of the bound-mode frequency reads 



IU>I 



Si; 
M 



Mfj, 



j + n+l 



(D13) 



To understand how this scales with M/x in the small-mass 
limit, we note that for a <C Mjj, and at first-order in P 



we have T[j + iP - a + 1] ~ -iP/P 2 
we get 



4M 2 fj, 2 



Finally 



Mut 



1 



iJ+l-s' 



(am-2r + n){MLiy J+ - i x 



tfl-^Tffi 



2]r[j 



(j + 1 + n) 2 ^+ 4 r[l + 2j] 2 r[2 + 2j] 2 T[n + l]T[-j + s'] ' 

(D14) 

The fundamental mode, n = 0, for the axial dipole (j = 
2 ,s' = 3) reads 



Mwr 



(o - 2r + /i) 



40(Af/i) 
19683 



li 



(D15) 



The formula above is valid when 0^ii< M/j, whereas, 
in the nonrotating case, it must be multiplied by a factor 
2 as explained above. A comparison with the numerical 
results for the non-rotating case and for the rotating case 
is shown in Fig. 0] and in Fig. [HI respectively. 
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We note the importance of the factor T[j + iP — a + 1] , 
which takes the form T[j — s' + 1] at lowest order in P 
and a/e and diverges because s' > j (s' = 3,j = 2). This 
is not the case in the axial perturbations of the Proca 
field (s' = l,j > 1) and the perturbations of the scalar 
field (s' = 0, j > 0) [iSSIialiil- It is this factor that 
contributes with a term (M n)~ 2S for the imaginary part 
of the quasibound frequency, resulting in a power-law of 
the form wi/n oc -(A/^) 4 ^ 2S + 5 = -(M^) 4 '+ 2S + 5 . 



^3. 

3 



-4 
-5 
-6 
-7 
-8 

-9? 

-10 



Proca (l=0,S=1,n=0) 
Num 

Ana 



0.06 0.08 0.10 0.12 0.14 0.16 0.1! 



1. Note on the monopole of Proca and massive 
spin-2 field 



by 



The monopole equation for the Proca field [79| is given 



d 2 u 



(2) 



dri 



2 , 2 2 6A/\ 



«(a; 



= 0. 



(D16) 
This can be written in the form (|D2[) taking a = 0, j = 
l + S = 1, and s' = 2. We can then solve analytically this 
equation in the same way as we did for the axial dipole 
and all the formulas apply. We then find that for this 
mode 



wj ^ 8(M^) 7 (n + l)(n + 3) 
T ~ (n + 2) 5 



(D17) 



in agreement with the numerical results of Rosa and 
Dolan [79J. In Fig. [8] we compare the numerical results 
and the analytical formula in the small-mass limit. 



FIG. 8. Comparison between the numerical and analytical 
results for the Proca field mode I — 0, n = as a function 
of the mass coupling M\i. The solid line shows the numerical 
data and the dashed shows the analytical formula ui/fi ~ 
-f(M M ) 7 . 



Another interesting behavior that we can infer com- 
paring with the axial dipole is that it seems that s' is 
simply given by the sum of the spin projection S and the 
spin of the field, i.e., s' = s + S = 1 + 1 = 2 for the 
monopole of the Proca field and s' = s + S = 2 + 1 = 3 
for the massive spin-2 field. 

Unfortunately the monopole equation for the massive 
spin-2 field (|30[) does not have a simple and understand- 
able form in the limit z <C 1 due to the complex form of 
the potential. However in the limit z>lwe can deduce 
the equation 



^2 ( a; - R o) 



;/ 
x 



6 



xR Q = . 



(D18) 



where Rq = <fio/ r - This looks exactly like the axial dipole 
equation in the same limit (|D3I) . By comparison we can 
see that the monopole acquires a centrifugal term with 
j = I + S = 2 in agreement with our numerical results. 
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